Hydrodynamics of fluid-solid coexistence in dense shear granular flow 
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We consider dense rapid shear flow of inelastically colliding hard disks. Navier-Stokes granular 
hydrodynamics is applied accounting for the recent finding that shear viscosity diverges at a lower 
density than the rest of constitutive relations. New interpolation formulas for constitutive relations 
between dilute and dense cases are proposed and justified in molecular dynamics (MD) simulations. 
A linear stability analysis of the uniform shear fiow is performed and the full phase diagram is 
presented. It is shown that when the inelasticity of particle collision becomes large enough, the 
uniform sheared fiow gives way to a two-phase fiow, where a dense "solid-like" striped cluster 
is surrounded by two fiuid layers. The results of the analysis are verified in event-driven MD 
simulations, and a good agreement is observed. 
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I. INTRODUCTION 



Granular shear flows have attracted much attention in 
recent years [ij. However, theoretical description of dense 
granular flows still remains a challenge The present 
work focuses on a rapid dense flow of inelastic hard disks, 
which undergo a shear motion. Here, the medium is flu- 
idized, particle collisions are binary and instantaneous. 
This brings about a possibility of a hydrodynamic de- 
scription of granular media ^ . 

There are still significant difficulties in hydrodynamic 
description of dense shear flows, even within the model of 
inelastic hard spheres. Although for low and moderate 
densities, a hydrodynamic description with the proper 
constitutive relations can be derived from kinetic theo- 
ries , constitutive relations for dense flows are presently 
unknown. Several attempts to extend the theory to high 
densities have been done recently 0, Jl]- It was 

proposed to apply free volume arguments [a, in the 
vicinity of close packing density, rimax = 2/(V3rf^), and 
use interpolation functions for constitutive relations be- 
tween low and high density cases @, iSi- At high 
densities, there is another problem related to the behav- 
ior of shear viscosity [1, [T^, [llj. While inelastic heat 
losses, thermal conductivity, and pressure diverge at the 
close packing density Umax, it was shown recently that 
shear viscosity diverges at a lower density Q . This may 
result in coexistence of fluid and solid phases [llj . 

We propose here novel interpolation formulas for con- 
stitutive relations between low-density and high-density 
cases, which account for viscosity divergence, and jus- 
tify them in MD simulations. We employ the resulting 
granular hydrodynamics for quantitative description of 
fluid-solid coexistence in shear granular flow and verify 
the hydrodynamic predictions in MD simulations. 



II. THE MODEL AND GOVERNING 
HYDRODYNAMIC EQUATIONS 

Consider a plane Couette geometry: an ensemble of 
inelastically colliding hard disks with unit masses and di- 
ameter d is driven, at zero gravity, by two walls. The top 
wall, located at y = H/2, moves in the a;-direction with 
velocity uq, the bottom wall moves in the opposite di- 
rection with the same velocity. The only parameter that 
characterizes the inelasticity of collisions is the coefficient 
of normal restitution, r. Boundary conditions in the x- 
direction are periodic. In the y-direction, no-flux and 
no-slip boundary conditions are implemented. Upon a 
collision with a driving wall, the normal particle velocity 
switches sign, while the new tangential velocity compo- 
nent is taken from Maxwell-Boltzmann distribution, with 
a mean equal to the wall velocity, and a variance corre- 
sponding to an instant temperature of the layer next to 
the wall. 

Now we present a hydrodynamic approach, which deals 
with the number density of grains, n(r,t), the granular 
temperature T{r,t), and the mean flow velocity [J]: 

dn/dt -f n V • V 0, 
n (dv/dt) = V • P , 

n(dT/dt) = -V • Q + P : Vv -r, (1) 

where P is the stress tensor, Q is the heat flux, and 
r is the energy losses due to the inelasticity of parti- 
cle collisions. The stress tensor P is given by P = 
[-p(n,r) -)-/i(n,T)tr(D)]I 4- 2r/(n,r)D, where D = 
(1/2) [Vt; -|- (Vu)-'"] is the rate of deformation tensor, 
D = D — i tr (D ) I is the deviatoric part of D, and I 
is the identity tensor, rj{n,T) and ii{n,T) are the shear 
(first) and bulk (second) viscosities. The heat flux Q is 
given by Q = —K{n,T)\IT, where K,{n,T) is the coeffi- 
cient of thermal conductivity. An additional term in the 
expression for heat flux, proportional to the density gra- 
dient [12, 13], can be neglected in the nearly elastic limit 
1 — <C 1 we consider throughout this paper. For higher 
values of inelasticity one also needs to take into account 
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inelastic corrections to transport coefficients [l2!|. How- 
ever, the main point of the present paper is addressing 
the challenging problem of describing a very dense flow 
and (as the first step) it is sufficient to assume the limit 
of nearly elastic particle collisions. 



III. TESTING CONSTITUTIVE RELATIONS IN 
MD SIMULATIONS 

For dilute and moderately dense granular flows, the 
Enskog-like constitutive relations are derived from kinetic 
theory The shear viscosity, r]E{n,T), the thermal 
conductivity KE{n,T), the inelastic heat losses TE{n,T), 
and the equation of state pE are given by 
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where Ge = 1^(1 - 7i//16)/(l - v)"^ and v = n{Trd^/4) 
is the solid fraction, and index E stays for " Enskog" . 
For small and moderate densities and for small inelastic- 
ities, transport coefficients are in a reasonable a gree ment 
with the results of event-driven MD simulations [ij] . For 
very large densities, one can use free volume arguments 

Q and show that all the constitutive relations (except 
for shear viscosity) diverge at the density of close pack- 
ing. However, there is strong evidence that at high den- 
sities the coefficient of shear viscosity behaves differently 
than other constitutive relations [3, diverging at a 
lower density than other transport coefficients [8|. An 
immediate consequence of this finding is a possible exis- 
tence of a solid-like phase, which is at rest or moves as a 
whole, as its density is higher than the density of viscos- 
ity divergence . Khain and Meerson fll| considered a 
very dense three-dimensional system and assumed a lead- 
ing order expansion of constitutive relations (except for 
shear viscosity) near the close packing density. However, 
to employ hydrodynamic equations within a wide range 
of densities, one needs a pragmatic approach involving 
interpolation of constitutive relations between the dilute 
and dense cases [1, S S HI ■ We will adapt this approach, 
suggesting several new interpolation functions, and test- 
ing them in MD simulations. 

To begin with, a global equation of state of hard disk 
fluid was recently proposed p = n T(l-|-2 G), where 

G = GE+m{i^rnax/ivmax-'^)-GE)- Here, mis an inter- 
polation function, given by m = [l + cxp((zyc — i')/too)]^^ 
with i>c = 0.70, nio = 0.0111, and lymax = 7r/(2-\/3), in a 
good agreement with MD simulations According to 
the same lines, we propose the modified form of inelas- 
tic heat losses, changing Ge to G in the corresponding 



expression in Eqs. ([2]). To test this expression in MD 
simulations, we consider a completely different system: 
a homogeneous freely cooling ensemble of inelastic hard 
disks in a box. In this case the third of Eqs. ^ becomes 
ndT/dt = — r, which gives the well-known Haff's law 
T = (1 + t/to)^^- Here temperature is measured 
in units of the initial temperature Tq, and time in units 

— 1/2 

of H Tq , where H is the system width. In our case, 
to = (7ri/V4) {d/H) [{l~r)G]~\ We measured the tem- 
perature as a function of time in MD simulations within 
a wide range of densities, while other parameters were 
kept constant. To ensure the homogeneous cooling state 
being stable [lB|, we took a sufficiently large restitution 
coefficient. As expected, we found the Haff's law to be 
vaHd and computed the value of to. Figure [T^ shows to 
versus f =< v > jvmax-, both from the theory (solid 
line) and MD simulations (circles). A good agreement is 
observed, in contrast to the Enskog predictions (dashed 
line), which deviate from the results of MD simulations 
at high densities. 



The next step is specifying the coefficient of shear 
viscosity. Consider a dense uniform shear flow. Here, 
the energy balance equation reduces to rj[du/dy)'^ = F. 
Measuring the temperature of the system in MD sim- 
ulations at different densities, we calculate the inelas- 
tic heat losses, and compute the shear viscosity rj. In 
Ref. Q, the coefficient of shear viscosity of a system 
of elastic hard disks was measured using the Helfand- 
Einstein expressions. It was found that rj diverges like 
a,, [vri — f^)^^7 where the viscosity divergence density is 
j^,, — 0.71, and a,, = 0.037 3, in contrast to the En- 
skog predictions. Our MD simulations confirm this re- 
sult. Figure [T}d shows that the Enskog formula (dash- 
dotted line) dramatically disagrees with MD simulations 
(circles) at high densities. An interpolation between di- 
lute and dense regime was also proposed in Ref. [1]: 
Vl = 77b [1 + o.n/{i'n — J^)], see Fig. [T}d, the dashed line. 
Since this interpolation is not accurate for the intermedi- 
ate densities, we propose an improved interpolation func- 
tion 7? = 77b [1 + a-qiy / v^Y / {vri — v) — a^/t',,]. The new 
formula (solid line) agrees well with MD simulations in a 
wide range of densities. 



In Ref. [3], the coefficient of thermal conductivity k 
was also measured in MD simulations. It was found to be 
larger than ke for intermediate densities {v ~ 0.55), and 
smaller than ke for larger densities, v ~ 0.75 (see Fig. 7 
in Ref. [81]). Thermal conductivity is known to diverge at 
the close packing density, as {vmax — . Incorporat- 
ing these findings, we propose the following interpolation 
formula k = ke{1 + Q-lv - lQv^° + Q.\l{v„iax - - 
.11 / Vmax) , which is in a good agreement with the re- 
sults of Ref. [1]. Now, we apply the resulting granular 
hydrodynamics to the problem of dense shear flow. 
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FIG. 1: Testing constitutive relations in MD simulations: 
heat losses (a) and shear viscosity (b). (a): A time scale for 
a temperature decay, to, as a function of density both from 
theory (solid line) and MD simulations (circles). The Enskog 
predictions (dashed line) do not agree with the results of MD 
simulations at high densities. The restitution coefficient is 
r = 0.999 (the simulation with the highest density was per- 
formed for r — 0.9999 and then the value of to was properly 
recalculated), the total number of particles is A'^ = 6480, and 
the system aspect ratio is 5 = 0.8. (b): The coefficient of 
shear viscosity normalized by the Enskog value as a function 
of density. Enskog formula (dash-dotted line) disagrees com- 
pletely with MD simulations (circles) at high densities. Also 
shown are the expression in Ref. (81] (dashed line) and the new 
interpolation formula (solid line). The total number of parti- 
cles is N ~ 6400, and the restitution coefficient is very close 
to unity (ranging from r = 0.996 to r = 0.99995 in different 
MD simulations) , to ensure that the solution of uniform shear 
flow is realized. 



IV. STEADY DENSE SHEAR FLOW: MD 
SIMULATIONS VERSUS HYDRODYNAMICS 



Let us measure the coordinate y in units of the system 
height H, the horizontal velocity u in units of the wall 
velocity uo, the temperature T in units of Uq, the density 



in units of 



the pressure P in units of rt 



7nax LiQ" 



Then the steady shear flow {d/dt = 0, w = 0, d/dx = 0) 
is described by the following system of equations: 



, /4T = const , 



ay V dy 



where R = (IB/tt) (1 — r) (H/d)"^ is the heat loss parame- 
ter, and the functions fi are the density-dependent parts 
of the constitutive relations: /i = k (-y/7rd/2)T~^/^, /2 = 
r]{2d^)T-^/^, lyG, and fi = n{l + 2G). No-flux 

and no-slip boundary conditions are given by dT/dy[y — 
-1/2) - dT/dy{y - 1/2) - 0, u{y - -1/2) - -1, 
u{y = 1/2) = 1. The total number of particles is con- 
served, which yields a normalization condition for the 
1/2 

density: /_j/2 ""(y) '^U ^ /' where / =< > / fmax is the 
average area fraction. Equations ([3]) allow for different 
solutions. The simplest one is the uniform shear flow. In 
this case, density and temperature are constant, and the 
velocity proflle is linear. For low and moderate densities, 
it is known that uniform shear flow becomes unstable, 
when the inelasticity of particle collisions exceeds a crit- 
ical value [1^, [13, [11] ■ This problem has been analyzed 
by linear stability analysis of the corresponding equations 
of granular hydrodynamics [H, [13, E| ■ The instability 
may result in shear-band formation, which was recently 
observed in MD simulations lip and analyzed theoreti- 
cally for moderate densities [I4I, see also 18 1. We found 
that for high densities, the same instability exists, but 
its threshold changes qualitatively. Indeed, as a result 



of shear viscosity divergence at z^^ < 



the uniform 



shear flow solution is impossible at sufRciently large av- 
erage density, / > v^flvmax- When this solution does not 
exist or is unstable, the density profile is no more homo- 
geneous. In this case, the regions with the density larger 
than that of viscosity divergence may appear. Although, 
those regions are at rest or move as a whole, the gran- 
ulate is fluidized there and granular temperature is not 
zero [Tl| . 

Our MD simulations show that when K is large enough, 
the system consists of three layers: an inner solid-like 
layer and two outer fluid layers, see example of a snap- 
shot in Fig. [2^. This fluid-solid coexistence has been re- 
cently observed in MD simulations [H, [l^l • To describe 
it hydrodynamically, one needs to solve Eqs. ^ for each 
layer separately, demanding continuity of the density, of 
the heat flux, and of the velocity at the interfaces between 
the layers. For the solid-like layer, the flrst equation of 
Eqs. ([3]) is replaced by u = const, whose value should 
be found from the overall solution. We find the solution 
of Eqs. ^ using a numerical shooting procedure, simi- 
lar to that described in [ll[. Interestingly, there exists 
a family of solutions (missed in (llj). which can be pa- 
rameterized by the position of density maximum in the 
system. In this study we only consider solutions with a 
symmetric density proflle (however, other solutions were 
found in MD simulations [l^ for moderately dense sys- 
tems). Consider a solid layer located in the middle of 
the system, so the density maximum is at ?/ = (as in 
Fig-Uti-)- FigureOa shows the density, temperature, and 
velocity profiles derived both from solution of Eqs. ^ 
(solid lines) and from MD simulations (circles). The den- 
sity (the upper panel) in the central part of the system 
(inside the cluster) is larger than the density of viscosity 
divergence. This results in m = plateau in the middle 
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(the lower panel). An excellent agreement between the 
hydrodynamic theory and MD simulations can be seen. 




FIG. 2: A two phase solution, (a): A snapshot of the system 
from MD simulations at t « 9 x 10"^. The solid- like layer in 
the middle (which is at rest) is surrounded by two fluid lay- 
ers, (b): The corresponding hydrodynamic profiles: theory 
versus MD simulations. The density (upper panel), temper- 
ature (middle panel), and velocity (lower panel) profiles in 
the system, both from solution of Eqs. (|3} (solid lines) and 
from MD simulations (circles). The total number of particles 
is A = 6480, the restitution coefficient is r = 0.99, the sys- 
tem dimensions are L ~ 97d, H ~ 80d. The corresponding 
hydrodynamical parameters are R — 323.86 and / — 0.7242. 



V. PHASE DIAGRAM 

Now we consider {R, < v >) phase diagram, see Fig.[3l 
Let the average density in the system be smaller than the 
density of viscosity divergence rt,, . If the hydrodynamic 
heat loss parameter R exceeds the threshold value (see 
the dash-dotted line in Fig. [3]) 



Rc 



{dfi/dv) hidfi/dv) 



2/3 



(4) 



the uniform shear flow becomes unstable (details of lin- 
ear stability analysis will be given elsewhere [2l[ ) . Above 
this threshold, different solutions with nonuniform den- 
sity and temperature profiles are realized. However, the 
two-phase solution is possible only when R is larger than 
some critical value R*{f) (computed from Eqs. ([3]), the 
solid line in Fig. [3]), so that the density contrast in the 
system is sufficiently large, and the maximal density is 
larger than n^. An example of this solution is shown 



in Fig. [2] (which corresponds to the upper asterisk in 
Fig. E]). In order to verify the hydrodynamic predictions, 
we performed MD simulations in different regions of the 
phase diagram. We found a uniform shear flow solu- 
tion (squares) below the dash-dotted line, while above 
the solid line, a two-phase solution is realized (asterisks). 
In the intermediate region (between the two lines) there 
are one-phase solutions with nonuniform density (and 
temperature) proflles. There are two types of symmet- 
ric (with respect to a density profile) solutions, which 
are allowed by Eqs. ([3]): the density maximum can be 
achieved either in the middle {y = 0), or near the walls 
(at y = ±1/2). MD simulations show that the system is 
denser in the middle when the average area fraction is not 
large enough (rhombuses), while for larger area fractions 
the density is maximum near the walls (circles). 
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FIG. 3: {R, < u >) phase diagram. The uniform shear 
flow is stable below the dash-dotted line, when R < Rc{f) 
[see Eq. @]. The two-phase solution is possible above the 
solid line, when R > R*{f) [where R*{f) is computed from 
Eqs. (|3}]. Symbols denote different solutions found in MD 
simulations: uniform shear flow (squares), two-phase solution 
(asterisks), nonuniform one-phase solutions with the density 
maximum in the middle of the system (rhombuses) or near 
the walls (circles) . The number of particles in MD simulations 
was N = 6450 ±30, and the system height was H = {79± l)d, 
while the system width L was varied in order to change the 
average density in the system. 



Interestingly, the solution with solid layer in the middle 
of the system seems to be metastable (see also Ref. To*] ) . 
After a sufficiently long time, the cluster starts slowly 
moving toward one of the walls until a steady state is 
reached. We followed the cluster dynamics by looking on 
the y-component of center of mass of the system, and on 
the velocity in the middle y = 0. Both quantities reach a 
plateau at long times, which corresponds to the position 
of the cluster near one of the walls. In this case, the 
cluster moves as a whole with a velocity which is slightly 
lower than the wall velocity. 
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VI. SUMMARY AND DISCUSSION 

The dynamics of a dense rapid shear granular flow in 
two dimensions is analyzed both theoretically (applying 
granular hydrodynamics) and by means of event-driven 
molecular dynamics simulations. We presented a general 
phase diagram and described different steady-state so- 
lutions. The most interesting solution describes a flow 
consisting of three layers: an inner solid- like layer and 
two outer fluid layers. This solution is possible due to 
the fact that viscosity diverges at a lower density than 
other constitutive relations [111, . A possible direc- 
tion of future research is increasing the streamwise length 



of the system. In this case, the layered structures may 
give way to wavy patterns [20]. For moderate densities, 
these structures can be explained theoretically [13] , while 
for higher densities viscosity divergence should be taken 
into account. 
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